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FOREWORD 


The Systems Technology Laboratory (STL) is a computational re- 
search facility located at the Goddard Space Plight Center of the 
National Aeronautics and Space Administration (NASA/GSFC) . The 
ST’j was established in 1978 to conduct research in the area of 
flight dynamics systems development. The laboratory consists of a 
VAX-11/780 and a PDP-11/70 computer system, alo.ng with an image- 
processing device and some microprocessors. The operation of the 
Laboratory is managed by NASA/GSFC (Systems Development and Anal- 
ysis Branch) and is supported by SYSTEX, Inc., Computer Sciences 
Corporation, and General Software Corporation, 

The main goal of the STL is to investigate all aspects of systems 
development of flight dynamics systems (software, firmware, and 
hardware) , with the intent of achieving system reliability whxle 
reducing total system costs. The flight dynamics systems include 
the following: (1) attitude determination and control, (2) orbit 

determination and control, (3) mission analysis, (4) software en- 
gineering, and (5) systems engineering. The activities, findings, 
and recommendations of the STL are recorded in the Systems Tech- 
nology Laboratory Series, a continuing series of reports that in- 
cludes this document. A version of this document was also issued 
as Computer Sciences Corporation document CSC/SD-81/6028. 

The primary contributor to this document was 

Joan B. Dunham (Computer Sciences Corporation) 

Other contributors include 

Anne C. Long (Computer Sciences Corporation) 

William Wooden (Goddard Space Flight Center) 

Single copies of this document can be obtained by writing to 

Keiji Tasaki 
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ABSTRACT 


This document, which Is an update to Computer Sciences Corporation document 
number CSC/SD-78/6002, describes the mathematical theory of the computa- 
tional algorithms employed in the Onboard Navigation Package (ONPAC) System. 
This system, which simulates an onboard navigation processor, was developed 
to aid in the design and evaluation of onboard navigation software. The mathe- 
matical formulations presented include the factorized UDU form of the extended 
Kalman filter, the equations of motion of the user satellite, the user clock equa- 
tions, the observation equations and their partial derivatives, the coordinate 
transformations, and the matrix decomposition algorithms. 

Use of the ONPAC system, with sample input and output, is described in a com- 
panion document, the ONPAC User’s Guide (Reference 1). 
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SECTION 1 - INTRODUCTION 


The Onboard Navigation Package (ONPAC) Simulator simulates an onboard nav- 
igation processor assembly using a modified version of the design proposed In 
Reference 2. The pseudorange and delta pseudorange observations for the es- 
timation of user spacecraft position and time will be measured by tlie onboard 
receiver/processor assembly from Information broadcast by the NAVSTAR/ 
Global Positioning System (GPS), The pseudoirange observation is modeled as 
the line of sight distance from the GPS to the user, and the delta pseudorange 
is modeled as the change in the pseudorange over a period of time. The systein 
is currently being used for analysis and evaluation of the algoi'ithms presented 
in this document for premission studies. 

The ONPAC system Is being developed on a Digital Equipment Corporation 
(DEC) PDP 11/70 computer which is similar to the DEC LSI 11 which the nav- 
igation processor assembly will use. One proposed use of the processor as- 
sembly is to be a part of an experimental package to be placed on Landsat-D 
for onboard orbit determination using Phase I GPS. A sample case and sample 
output for Landsat-D are presented in the appendixes of the user's guide (Ref- 
erence 1). 

1.1 OVERVIEW OF ONPAC CAPABILITIES 

The ONPAC system processes GPS pseudorange and delta pseudorange obser- 
vations sequentially to estimate and apply corrections to a host vehicle state, 
which Includes the satellite position and velocity, two terms to describe cor- 
rections to the host vehicle clock, and a satellite drag coefficient. The inte- 
grator used to predict the position and velocity from observation to obseiwation 
is an Euler Integrator; the force model used In this Integrator may be varied 
by the user. The force model options are detailed in Section 3.2. The covari- 
ance matrix Is propagated with a state transition matrix which is a Taylor 

9 3 

series expansion of the analytical state transition matrix to At , At“ , or At 
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as specified by the user. A state process noise covariance matrix is also com- 
puted to allow for errors in the knowledge of the state and is added to the co- 
variance matrix. The variance of the state noise in position, in the drag 
coefficient, and in the clock terms may be specified by the user. 

The ONPAC system design is modular so tliat algorithms can be replaced or 
added with minimal changes to the existing simulator. Procedures for maldng 
such changes are discussed in the user's gttide (Hef 'rence 1). 

1.2 SUMMARY 

This document describes the basic mathematical algorithms used in ONPAC. 

Derivations for most of the algorithms are available in standard texts and 

sources are given in the references. Section 2 describes the extended Kalman 

T 

filter and gives a brief derivation of tlie UDU filter. The step-by-step appli- 
cation of the UDU*^ filter in ONPAC Is presented, with the points identified at 
which editing and smoothing are performed. The state vector and the state 
propagation equations are given in Section 3. The mathematical models in 
ONPAC use four different coordinate systems. Transformations from the 
inertial to the other three systems are also given in Section 3. The pseudo- 
range and delta pseudorange observation models and partial derivatives are 
described in Section 4. The state process noise covariance matrix used in 
ONPAC is given in Section 5. Appendix A describes the matrix decomposition 
algorithms. The data simulation which is performed using the Goddard Tra- 
jectory Determination System (GTDS) is described in Appendix B. The relation- 
ship of these elements of the ONPAC program is shown in Figure 1-1, 
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SECTION 2 - ESTIMATION EQUATIONS 


Matrix factorization teclmlquos have been demonstrated to improve the stability 
and accuracy of Kalman filters (Reference 3), The factoring of the covariance 
matrix, P, into square root 


T 

P « SS 

or upper trianguiar, U, and diagonal, D, components 

P = UDu'’' 


reduces the occurrence of numerical problems by keeping the covariance ma- 
trix positive definite. The filter in ONPAC is an application of the U-D for- 
mulation of an extended Kalman filter (EKF). 

2.1 EXTENDED KALMAN FILTER 

A derivation of the EKF, also lmo\vn as the extended sequential filter, can be 

found in many sources (Reference 4 and 5).. The EKF equations are presented 

T 

in this section as a basis for the UDU filter discussion which follows. 

Given a state vector and covariance matrix at time t^^ ^ , the prediction 
and update equations at time are as follows; 

1. Given X(t, J , the estimated n-dimensional state vector at t, , , 
^ - k-1' k-1 

and covariance at t^^_^ , integrate the state equations of motion 

t), 


X = F (X, t) 


(2-1) 
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from to , with the Initial conditions » This gives the pre- 

dicted state X(tj^^ , The carets over X and P Indicate the updated estimates 
after processing an observation. 

Compute the state transition matrix 'if (t^^, , by integrating the differen- 

tial equations 

tj^_i)-A(t)$(t, (2-2) 

where A(t) Is the matrix of partial derivatives 


A(t) ^ 


'ap(x, t) 
dx 


(2-3) 


evaluated at x “ Initial conditions for the integration are 

Propagate the covariance matrix using the state transition matrix to obtain the 
predicted covariance matrix at time tj^ (without process noise), 

2. Obtain tlie observation at time tj^, T(tj^ • Compute the predicted 
observation, G(X(tjp, tp , where G is a nonlinear function of the state pa- 
rameters and time. 

Compute the (1 X n) matrix of obseivation partial derivatives 


tj^_^) - 1 , the Identity matrix. 


H(tk> = 




J X=X(t, ) 
“ K 


(2-5) 
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The observation residual, or observed minus computed n??servatlon (0 - C) Is 
given by 


g(V - Y(tj^) - G(X(y . y 


( 2 - 6 ) 


3. The Kalman gain, K(tjp , Is 


K(y = p'^(y n(t,^) p(y 


(2-7) 


2 

whore <j Is the measurement variance, 
m 

The quantity within the brackets Is a scalar and K(y is a vector with n com- 
ponents. 

The Kalman gain is usev’ to compute the updated state vector 


X(y -X(y +K(y g(y 


and the updated covariance matrix 

p(V = [i-K(y H(y] p(y 


( 2 - 8 ) 


(2-9) 


The updated state vector and covariance matrix at time tj^ are then the input 
for propagation to time in the next step, 

2.2 UDU^ FILTER DERIVATION 

The covariance matrix is a positive definite square matrix which may be factored 
into a triangular matrix and its transpose. It may also be factored into a trian- 
gular matrix with imity on the diagonal, its transpose, and a diagonal matrix. 
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Tho EKI^ equattons can then be formulated using the component matrices instead 
of the covariance mattiX, Tho covariance matrix is computed as shown below. 

P UDu"^ (2-10) 


where U is an upper Ulangular matrix and D is a diagonal matrix. 
The Initial conditions for the tj^th step of the EKF filter are 


°<‘k-i> 

Propagating the covariance matrix (^vithout process noise) to time tj^ yields, 


P(tk) = D(y u (y 

= ‘k-i> ^<‘k-i> ®<‘k-i> ‘k-i> 


( 2 - 11 ) 


*<‘k- ‘k-i> °(‘k-i>®<‘k-i> ‘k-i> ®<‘k-i>] 


T 


The covariance update equation is 

P(y = [I-K(tj^) H(y ] P(y (2-12) 

For the following equations, all quantities are evaluated at . Replacing 
K in Equation (2-12) with Equation (2-7) gives 

? = P - P^h'^ ( HPh’^ + HP (2-13) 
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Substituting Equation (2-10) Into Equation (2-13) gives 


A a/s/nT 

P ° UDU 


VD\F - UDU V ^HUDU^H^ + HUDU'^ 

T T / T T 2 \ ^ T 

U D - DU H (HUDU’^H 1 HUD U"^ 


( 2 - 1 * 1 ) 


For a single observation, the term 


hudu’^h'^ + 


is the scalar 


1=1 \ j-i ' 


(2-15) 


where n is the dimension of the state vector, d. = d. . , the diagonal terms 

X ^9 ^ 

of D , h. = h. . , the elements of the U matrix and u. , are the elements 

11,1 1, j 

of the U matrix. 

Substituting this relationship into Equation (2-14) gives 




D 


(duVKdu^h^)' 


a 


n 


U 


(2-16) 


-1 -T 

Equation (2-16) is then premultiplied by U and postmultlplled by U to 
give 





(2-17) 
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L 


if 


II 

<j 


n 
; { 


vf 


</ 




(j 


jrst:? 



4 

!l 

4w.:;’ 


which Is, substituting in (2-16) 


= D - 


(du’^h’’^)(du’'^h'^) 


a 


(2-18) 


/S -1 -1/N -X 

Because both P and U (and U ) are positive definite matrices, U ‘PU is 

also positive definite (Reference 3). Therefore, the matrix (2-18) is a positive 

definite matrix and can Itself be decomposed into upper triangular and diagonal 

component matrices. 

Let 


M =D - 


T T / T T\ 
(DU H ) (dU H ) 

0!„ 


/o 


^19) 


and 


M = BWb’^ (2-20) 

where B Is an upper triangular matrix with unity on the diagonal and W is a 
diagonal matrix. 

Then, 


/N/N/nT T T 

UDU =UBWB U 


and the following can be identified; 


U = UB 


D = W 


(2-21) 
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rr 


II 

i 




I 


I 

I 

!t 


The matrices B and W can be found by applying the factorization equations 
given In AppendLs: A. 


The general terms for components of the matrix M are 



( 2 - 22 ) 


and for i 7^ j 



(2-23) 


: i 

I 

V £-* 




After decomposition and some manipulation, the general term for W is 


w. = 
1 


d. a. . 
1 i-1 

% 


where w, = w, . 

i 1, 1 


a 


i-1 





(2-24) 


(2-25) 


and 



(2-26) 
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The general off-diagonal term of the B matrix is 


>.(vi-|\v.-,)(vS\v,l 

i-j > 1 iv ~ ' 

“i-1 


(2-27) 


Then, the update equations for the U and D matrices r 


are 


“1 = 


d, ot 

a 


- i 1-1 


1 


and 


(2-28) 


d. S S, u d. S 

u, =U t J ® 

i.i i.j <T~ ,4-, — 

j-1 Is^l 


k J 


a 


3-1 


(2-29) 


where S Is defined as 




(2-30) 


The Kalman gain Is given by 


K = U DU*^ ( H UDU*^ 

\ m 

UDu"^ 


^-1 


a 


n 


(2-31) 
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The components of the gain vector are 



(2-32) 


2.3 ONPAC FILTER 

The filtering algorithm used to process each observation in ONPAC is discussed 
below. f 

1. Enter with a state vector, X(tj^ , and covariance matrix compo- 

nents, U(t^^_^) and D(tj^_^) , from the previous step. If this is the 
initial step in the filter, initialized. 

Sample initial conditions are given in the user's guide (Reference 1). 

2. Retrieve , (the time of the next observation), and Y(tj^ , (the 

observation). Correct the time t^ with the previous estimate of the 
user clock error (Section 3.3). Propagate the state to t^^ using the 
Euler integrator (Section 3). If At = t^ - ^ is larger than the 

maximum allowable stepsize, the integration from ^ to t^^ is 
done in substeps no larger than At 

max 
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3, ComputG (tl^6 state transition matrix) and Q(tjP (the 

noise matrix), 

4. Propagate the covariance matrix 


Ct.' ‘k-l* ‘k-l> ^ 


5, Decompose P(tj^) into U and D components using the method given 
in Appendix A 


. p<y = u(tj^) D(y 


(2-34) 


6. Fade the filter memory by multiplying the D component of the fac 
tored covariance matrix 


D(tj^) = p D(y 

where 

4 

p = i V i3n(i) + /3An(i) (/3 . <p</3 ) (2-35) 

i=l 

where and are vectors of smoothed residuals from each of 

the GPS satellites in the current constellation, and B . and B 

min '^max 

are the minimum and maximum values allowed for p . If there are 
fewer than four GPS satellites in view, this step is not performed. 

7. Compute the observation 


G(x„. y =R(v 


for pseudorange 


= AR(tj^) for delta pseudorange 


(2-36) 
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and the matrix of observation partial derivatives 





for range 
for delta range 


These equations will be given In Section 4. 

8. Compute the observation residual 


g<v = ’'(V - ■ V 


9. Compute the a's and S's 


OL 


i 



+ 



S. = h, + 
J j 



htf u 


i 


j 


(2-37) 


(2-38) 


(2-39) 


(2-40) 


where 


2 


m 


2 

= CTj:^ for a pseudorange observation 
= for a delta pseudorange observation 


10. Test the residual for acceptability by computing the square-to- 
variance ratio 



(2-41) 
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11. Test p~ against a tunable parameter. If pj^ ^/^^nax ’ 

the observation is rejected. In that case, set 

= X(y 

S(‘k-i) = 

‘k-1 " ‘k 

^ ” Wl 

and go back to st'ijp 1, 

If p < p and the observation is a range, compute the smoothed 
K nicix 

residual 

I3^{1) = +p(pj^- (3jj(i)) (2- 

t 

If the observation is a delta range, the smoothed residual is 

^AR^^^ + P (\ “ ) (2- 

where (/3j^(i), is the smoothed square-to-variance ratio from 

the previous observation of the ith GPS and p is a tunable parameter. 
After is computed, it is checked against a tunable 

parameter P£ and, if larger than P 2 , set equal to p£ • 

12. Update the U and D components of the covariance matrix and com- 
pute the Kalman gain, K, using Equations (2-28), (2-29), (2-30), 
and (2-32). 


42a) 


42b) 
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13, Update the state vector 


x(y -s(y •*’K(yg(y 


(2-43) 


If there are more observations, go to step 1. 
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SECTION 3 - STATE PARAJvIETER PROPAGATION EQUATIONS 


The ONPAC simulator has nine stato paramotors, the position and velocity 
components, two clock terms, and a drag coefficient. For satellites such as 
Landsat, which are not highly drag perturbed, this can be reduced to eight by 
leaving out the drag coefficient. 

With the Intention of reducing the computation time and storage needed, the 
Taylor series expansion in the algorithms for the stato pi'opagatlon were trun“ 
catod to the minimum number of terms necessary for achieving the desired 
accuracy. Observations will be made very frequently at the rate of a pseudo- 
range and delta pseudorange pair in every 6.6 seconds when possible. This 
will reduce the Impact of the neglected terms in Uie propagation algorithms on 
the filter accuracy. 

3.1 STATE VECTOR 

The state vector, X , which Is used in ONPAC is given by 


X 

y 

z 

b 

X 

y 

z 

b 

d 

•m m 


(3-1) 


whore (x, y, z, x, y, z) are the Cartesian position and velocity components 
in Earth-centered Earth-fixed (ECEF) rotating coordinates, (b , b) are the clock 
bias and bias rate expressed In kilometers (km) and kilometers per second 
(lon/sec) (l.e. , the bias and bias rates are multiplied by the speed of light 
yielding Ion and km/sec. This provides immediate comparison with the errors 
in pseudorange and delta pseudorange) and d is the drag coefficient 


d ~ 


C_A 
D cx 

2 m 


(3-2) 


where A ~ the cross-sectional area of the satellite 
cx 

m » the mass of the satellite 
= a constant coefficient 

3.3 FORCE MODEL 

The ONPAC force model is variable, to be set by the user at run time. The 
equations of motion include the force terms due to the central body attraction 
and the coriolis and centripetal terms. To this may be added the force terms 
due to J, or to a 2x2, 3X3, 4x4, or 5X5 geopotential model and 

A 

atmospheric drag. The gravity potential is a Pines model (Reference 6). 

The gravity potential function uses a four parameter model. The potential 
function Is 




nmax/R 

® 


,n 


E 


n=l 


J A , 
n n, 0 


(u) 



r 

m 


(s,t) +S 


n,m 


i 

m 



(3-3) 
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where 


C 

n,m 


R mean ratUus of the Earth 

0 

jLi o GM, the gfavltationi:^ parameter (the gravitational con- 
stant times the mass of the Earth) 

J c zonal harmonic coefficients 
n 

S = tesseral harmonic coefficients (the values used are from 
GTDS) 

9 o 2 l/2 

r £3 (v* + + z") £3 magnitude of the position vector 

s = x/r 


t " y/r 
u ~ z/r 

A (u) a d^ ^n(’^)/'^u^ (where Pj^(u) = the Legendre polynomial 
of the first kind of degree n) 


nmax “ a user specified input which =0, 1, 2, 3, 4, or 5 
= 0 to pr’oduce a 2 -body only geopotential 
= 1 to produce a Jt geopotential 
s 2, 3, 4, or 5 giv/ng a 2X2, 3X3, 4X4, or 
5X5 geopotential field, respectively 


The functions r^(s,t) and i^(s,t) , the real and complex portions of the 
potential expansion, are defined by 


(s + j t)^^ = r^(s, t) + j i^^(s, t) 


where j = , 


The gravitational acceleration is the gradient of the potential function with re- 
spect to the four parameters (r , s , t , u) 


3V 3v 3V oV 

F Ar+*^As At + ~Au 
3r 3s 3t ou 


(3-4) 
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The acceleration in Cartesian coordinates duo to the gravitational potential is 


a “ a, “ sa . 
X 1 4 


a « a„ - ta . 
y 2 4 


a s a„ - ua, 
z 3 4 


a, « 


1 &V 


1 r bs 


„ 1 sv 
°'2 " 7 bt 


■ 

^3 r b\x 




The equations of motion for the spacecraft position vector are 


(3-5) 


dx 


dt 

at 


= X 


(3-0) 


dt 


The equations of motion for the spacecraft velocity vector are 


+ 2fiy + - d?7vx 


‘Ht’ ” % ” 


(3-7) 
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dt 


aa 

z 


(.Irjvz 


whoi’o d tho drag seale faGtor 

2 

„ Cn X satellite area in km 
^ 2 X satellite mass in kg 

Cp “ coefficient of drag 

222 1/2 

V » (X + y + z ) = magnitude of velocity vector 

Sk ® Earth's rotation rate 

t\ S3 the atmospheric density, computed as follows; 

The height of the satellite, h , above the reference ellipsoid is found from 


h = r 



(3-8) 


where € =* Earth's elUpticlty 

S3 Semi-major axis of tho reference ellipsoid. 

The height Is scaled with 


S Tq + ll 


(3-9) 


where Is a constant. 

The scaled height h Is compared against a series of heights 

s 


/ill = h - h (k) , fork -1,6 (3-10) 

s \ 

where h. is a vector of threshold heights. 

V 
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Wlion the first value of k Is found for wlileli 


Ah>0 

the atmospheric density is computed as 

rj « W^(k) [w^(k)4h + x] 

where W , W, , and W arc vectors of constant‘s, 
a b e 

If 


then 




'1 = 0 . 

The equation to model the behavior of the drag coefficient is 



where is an input parameter, the time constant of the drag. 

3. 3 CLOCi^w MODELING 

The clock equations are as follows; 


db ,• 

5T='^ 

dS 


where is the time constant of the clock. 


(3-11) 


(3-12) 


(3-13) 


3-0 



QRSGtNAi- ’t 

OF POOR QUM-'TY 


For a stop of At from to tj^ , 


(3-14) 

»5(V = 5^(tk_i) 


3.3.1 Correction, of the User Clock 

The independent variable, t , is the user clock. The corrections to the user 
clock, (b , lo), which are estimated in the filter, must be applied to this in- 
dependent variable. 

When an observation is received, the time at which it is received must be 
corrected with the predicted offset and rate. The clock rats equation (Equa- 
tion (3-9)) computes the correction in kilometers. To convert it to seconds, 
it needs to be divided by the speed of light in kilometers per second (c). 

The true time, at the kth observation, is approximately 
\rue ^obs 

w t - 
obs 


;k 

c 


b 


k-1 u. ^ ^ ^k-1 

true k-1 c 


or 


k-1 


obs 


+ 1 


k-1 


k-1 


true 


k-1 


(3*T5) 
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Tho observed time Is corrected to be t^^ by 


t 


k-1 


t 


bbs p 


1 + 


+ t 


lc-1 


k-1 c 


~'c 


(3-lC) 


3.4 MODII’IED EULEIl INTEGIUTOR 

The slate equations for the spacecraft position and velocity are Integrated using 
a modified Euler Integrator. The procedure is described below. 

Using the state vector and acceleration vector from the previous step, at , 
compute the position and velocity at 







4- At 
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v“ S3 V 

~k -k-1 




a/ .. 

“ 2 -k-1 


(3-17) 


whore At is the integration stepslze. 

This now state i’p is used to compute the acceleration at tj^ , r^^ by 
calling the dox'ivative evaluation subroutines. Then, the two acceleration values 
(r, , r, ) are used to perform the integration 

'~lv-l — IC 


• • tilt ft. * • . 

r, = r, - -7T fe 1 ) 

-k -k-1 2 ^k-1 -k 


^ •• 

4v~-k-l ‘ “''-k-1 3 -k-l‘" 6 -k 


r, = r, + At r. 


(3-18) 
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The stepslze, At , is moiiitorod so that 


At ^ At 


max 


whore At is an input to the Integration routines, 
max ^ ® 

3.5 COORDINATE TltiVNSFORMATIONS 

The GPS satellites broadcast their Inlormatlon In Earth-centered Earth-fixed 
(ECEE) Cartositm cocrdlnatos, tmd the ONPAC ephomerls and estimation al- 
gorithms are all written using ECEF coordinates. More information on tlie 
filter behavior and ori'or sources can be obtained from examining ONPAC re- 
sults in otlier coordinate frames (in particular, the UVW coordinates and the 
Keplerian orbital elements). To convert from ECEE coordinates to UVW coor- 
dinates or to Koplerian elements first requires a conversion to an inertial 
coordinate frame. The one used for this purpose in ONPAC is the Eartli- 
centered inertial (ECI) coordinate system. 

3.5.1 Eartli-Centered Inertial (ECI) to Earth-Centered Earth-Fixed (ECEF) 
Coordinate Transformation 

In general, the transformation from an ECEF to an ECI coordinate S 3 'stem is 
obtained by the rotation 


r 


"COS Ci t 
sin Q t 
- 0 


-sin 0 t 
cos a t 
0 


0 

0 

1 


r' 


(3-19) 


where r is tlie satellite position vector in true of date ECI coordinates, r' 
is the position vector in ECEF coordinates, Q is tl’.e EartlFs rotation rate, 
and t tlie elapsed time from the epoch of the tnie of date sj-stem. For ONPAC, 
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it is necessaiy only to Icnow the instantaneous correction from one system to 
the other. In tliis case, the time from epoch is zero and the transformation 
from ECEF to ECI is 


r = r‘ 


r = r' + O X r 


(3-20 


where Q> is the vector directed along the Earth's North polar axis whose mag 
nitude Is , the rotation rate of the Earth. 


3.5.2 ECI to UVW Coordinate Transformation 
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The transformation from ECI to UVW is 



/s 

u • i 

A A 

u • 3 
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u • k 
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/N 0* 
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V • 3 

A 

V • k 
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A A 

A 


w* i 

w* 3 

w • k 


(3-24) 


where (i', j", 1c) are the ECI unit vectors and r" is the position vector in UVW 
coordinates . 


3.5.3 Keplerian Orbital Elements 

The orbital elements are determined using ECI satellite position and velocity 
(r, r) . The angular momentum is 


e = r X f 


= c.i + c.j + c, k 
1 ] k 


The vector along the nodal line is 


n = k X c 


= n.i + n .3 

1 y 


O' , o 

= -c,l + 0.1 
I 1 


(3-26) 


(3-26) 


The Laplace vector, whose magnitude is the eccentricity, e , and which points 
along the periapsis line is 


e 



r - (r • r) r 


e = 



(3-27) 
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andjj, = GM , the gravitational parameter. 
The semilatus rectum is 


P == 


£ » S. 


(3-28) 


The semimajor axis is 


a = 


P 

1 - 


2 

e 


The inclination of the orbit to the X-Y plane, I , is given by 


(3-29) 


cos I = 



(3-30) 


The inclination is always less than 180 degrees. 

The longitude of the ascending node, measured from the X axis, is given by 


n. . 

cos Cl = 777 (3-31) 

ISl 

If m <0 , the orbit is retrograde and the longitude = O + 180° . 

The argument of perigee is given by 


. n • e 
cos 00 = 

n 1 e 


(3-32) 
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If < 0 » the argument of perigee = w + 180°, 

The true anomaly Is given by 


If r • r < 0 , the true anomaly » f + 180° . 

The eccentric anomaly E Is given by 

< 3 - 34 ) 

3. 6 STATE TRANSITION JMATRLX 

A state transition matrix propagates the estimated corrections to the state vsc= 
tor for\vard in time and it is computed by integTating the linearized form of the 
equations of motion, In ONPAC, tlie covariance mati’Lx is propagated by using 
an analytic approximation to tlie state transition matrix. The analytic state 
transition matrix is expanded in a Taylor series and may be truncated at 
0(t), 0(t ) or 0(t‘ ) at the user's option. The force model for the state trans- 
ition matrix includes the central body attraction and the corioUs terms. The 
user may also include some terms due to drag. 

The general differential equation for the state transition matrix $ is 

i(t, tj._j_) =A(t)$(t, tj^_^^) (3-35) 
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which has the solution 


(t 


k’ Sc-1^ 



A(t) dt 


« exp [A(t^__j^) At] 


(3-36) 


2 2 
A (t .) Ar 

+ ... 


where At = t^ - tj^_^ . The matrix A(t) Is composed of the partial derivatives 


of the equations of motion 

F(t) , with regard to the state, 

X . 



°3X1 

^^12 

^3x1 

’'*13 


°1X3 

1 

°1X3 

At 

0 

O 

1 

11 

hi 

^3X1 

^22 

°3xl 

^23 


°1X3 

0 

^1X3 

1 - At/ Tj 

0 


°l!<3 

0 

o 

f-* 

X 

03 

0 

1 - At/r^ 

where 0i^ = I (and is a 3 by 3 

submatrix) 




0 . = 3 by 1 zero submatrix 

0._ = lAt + 0 At"'/2 (and Is a 3 by 3 submatrix) 

Qi 

0 = 1 by 3 zero submatrix 

1 AO 

021 -0^ At (and is a 3 by 3 submatrix) 

0 = I + 0 At (and is a 3 by 3 subniatrix) 

AM 


At = 




(3-37) 
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(3-41) 


where O = rotation rate of the Earth. 

The equations of motion, F(t) , are Equations (3-6) and (3-7) defined in 
Section 3. 2. 


3.7 FADING MEMORY 


The most recently made observations are weighted more heavily than 
previous ones by fading the filter memory. When the covariance is pro- 
pagated, as described by Equation (2-4), an additional factor is included 


P(U = 

Iv 


t JP(t, 


ic k-1 k-l k k-1) 


T 


(3-42) 


3-15 



whoro the variable s can be described as 


s = 0 (3-43) 

with ^ the interval between observations and r a time constant. Since the 
covariance Is factored into U-D components for gain and update computation , 
the D component of the propagated covariance Is multiplied instead of the 
entire covai’Iance, 

I 

D(y - sDiy (3-44) 

which is equivalent to (3-42) when there is process noise covariance. In ONPAC, 
the s term is approximated with a combination of smoothed square-to-variance 
ratios of range and delta range observations from the GPS satellites in the 
constellations, as described in Section 2,3. 

The fading memory can be compared to the process noise covariance matrix 
given in Section 5. The process noise terms, which are added to the covariance, 
provide a minimum value for specific terms of the covariance at any given time. 
The fading memory multiplies tlie entire covariance and can impact the covari- 
ance much more than the process noise. This limits the use of the fading 
memory to periods of good GPS satellite visibility (four or more in view). 

When used together, the fading memory will overwhelm the proce >s noise but 
there is reason for having both available to ONPAC. When the GPS satellite 
visibilitj' is poor, only the process noise is used. When it is good, both the 
process noise and the fading memory are used. 
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SECTION 4 - OBSERVATION MODEL AND PARTIAL DERIVATrVES 

The pseudorange observation equation Is the equation for the length of the llne- 
of-slght vector from the user to the GPS satellite, The delta pseudorange ob- 
servation equation Is the difference between two range measurements made over 
a short time span. The partial derivatives of the pseudorange equation are 
straightforward. Those of the delta pseudorange are the first terms of an ex- 
pansion In which the assumption Is made tliat a linear approximation will suffice 
over a short time span. 

4.1 PSEUDORANGE OBSERVATION MODEL 

The pseudorange observation equation at time tj^ Is 







,11/2 


+ b 


(4-1) 


where (x, y, z) - the user satellite position at t^^ 

(s , s , s ) = the GPS satellite position at t, 

X y z' ^ k 

b = the user clock bias at t, 

k 

The partial derivatives of the pseudorange observation with respect to the state 
parameters are 


BR 
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oz 




db 


= 1 


(4-2) 
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(4-2) 

(Cont'd) 


where the components of the unit vector along the line-of-sight from the user 
satellite to the GPS satellite are given by 




s (t, ) - x(t, ) 
x' k' ' k' 


''x^V R(tj^) - b(y 


A„(t,) 


S - y 

y 


y'k' R-b 


(4-3) 


s - z 
z 


z' k' R-b 




4.2 DELTA PSEUDORANGE OBSERVATION MODEL 
The delta pseudorange observation equation at time tj^ is 


AR(y = R(y - R(tj^_^) 


(4-4) 


This equation assumes that a pseudorange measurement precedes the delta 
pseudorange measurement at time 

The partial derivatives of the delta pseudorange AR with respect to the state, 
X , at time t^^ are 


3AR(t|P 5R(y SR(t|^.j) 

sx(y ° ax(y ■ sx(t^) 3x(y 


3H(y 




(4-5) 
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The state transition matrix, can be approximated by the expansion 

For the nine-parameter state, a further approximation for the matrix A can bo 
introduced, where 
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Then, the partial derivatives with respect to the state parameters are 
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oAR 

ad 


0 


whore At *=» tj. - • In the pseudorango measurement, the At Is small 

(on the order of 0. 6 second). The ONPAC program assumes that a delta pseudo- 
range observation at time t. is following a range observation at t, . 



SECTION 5 ~ DERIVATION OF THE STATE PROCESS NOISE 
COVARIANCE MATRIX 


The covariance matrix in a Kalman filter or EKF will become saturated, the 
terms becoming very small, after a large number of observations have been 
processed. The result is such that the filter will no longer significantly correct 
the state (i.e, , the state corrections become infinitesimal). To prevent this, 
a noise matrix which compensates for neglected terms in the force model is 
computed and added to the covariance matrix when it is propagated. Then, 
Equation (2 -4) becomes 




(5-1) 


where Q(tj^) is the state process noise covariance matrix. 

The noise errors are assumed to be uncorrelated in time (see Reference 4) and 
the Tvosition-velocity and drag terms are uncorrelated with the clock terms. 

The vector X is defined as the error in the state. 

-e 


6b 


X 

~e 


6b 


L6dJ 


(5-2) 


where 6r = position error 
6b “ clock bias error 

6v = velocity error 

« 

6b = clock drift rate error 
6d = drag scale error 
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Then, the differential equations to describe tliis system can be written 




+ w (t) 


( 5 - 3 ) 


where w(t) is a random forcing function to account for the errors in the force 
models used in the differential equations. This is called the state noise. The 
function w(t) is assumed to satisfy the following; 


EU(t)] = 0 

E w(t) w^(s)j = W6(t - s) 


( 5 - 4 ) 


where 6 = the Dirac delta function 

6(t - s) = 1, t = s 
0, t 7 ^ s 

E = the expectation operacor 
W = spectral density matrix 

This model assumes that the state noise Is unbiased and uncorrelated in time. 
The covariance matrix for the noise is 


Q(t) = E 




J ^ J ^ U) w(u) (s) (tj^ s) du ds 


k-1 k-1 
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= j f <l>(tj^, u) E w(u)u;'^ (s)J «l'^(tj^,s) du ds 
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If tlio position-velocity and drag terms are uncorrolated ^vith the clock terms, 
the state noise matrix can be partitioned and the clock terms derived separately. 
The noise matrix can be divided into the following: 


Q(t) 
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The clock process noise covariance matrix terms are 


q, 


= |c?At + 2(T^At^ + |7r2c?At^ 


bb 2 0 
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where 


and 


I’d is 

o 

r“ is 

0 . 
r; IS 


r2 is 


T 


min 


min\ min/ 

the clock white noise Allan variance factor 

the clock flicker noise Allan variance factor 

the clock random walk Allan variance factor 

the deterministic clock drift rate variance 

is the scaled minimum time where flicker noise predominates 
on the Allan variance ciirv’e 
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Q'(t) is defined as the subset of Q(t) which concerns Uie position-velocity 
and drag terms 
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Equation (5-4) can be rewritten In the form 
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6(u - s)J 


2 2 

for the position- velocity and drag noise terms where and cr^ are the rates 
of the noise variances for the acceleration and the drag coefficient. 

T 

Because Equation (5-5) for computing Q(t) contains bs^th and ^ , the state 
transition matrix, <l> may be approximated with fewer terms than those given 
in Section 3. 

The transition matrix is approximated for use in Equation (5-5) by 
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where At = t, - 1, , 
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Q = rotation rate of the Earth 
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(5-16) 


2 2 2 3 

The terms cr^ At and At were neglected since they are smaller than the 
neglected terms in the state transition matrix. 

The state noise process covariance matrix is a positive definite matrix. Con- 
sideration of the position-velocity components CQ--> Qir>3 shows that 

X X Xm 
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is a matrix which may be factored into U - D components. Let Q" be the subset 
of Q described above. Then, 
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whore 
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and q is a diagonal 6 by 6 matrix with elements 

A 
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APPENDIX A - MATRIX DECOMPOSITION 
A. 1 CHOLESICSf DECOMPOSITION 

Any positive definite square matrix can be factored Into a triangular matrix 
and Its transpose. If A Is such a matrix, 

T 

A = BB 

The matrix B Is not unique and may be either an upper or lower triangular 
matrix. If It Is an upper triangular matrix, the decomposition algorithms for 
an n X n matrbc are as follows: 
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for = j - 1 to 1 
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A. 2 UDU^ DECOMPOSITION 

The matrLx A can also be factored into a triangular matrix with unity on the 
diagonal and a diagonal matrix, 




A = CEc'^ 
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which can also be written as 


A = (CE^^^)(CE^^^) 


Implying that 


B = CE 


1/2 


Wlien C Is an upper triangular matrix, the algorithms for the C and 
ponents are as follows; 
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for i != 1 to n - 1 


Then, for j - 1 to n - 1 the diagonal terms are 
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and the off-diagonal terms in C are for m = 1 to n - j - 1 
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The diagonal terms of C are 


c. .-I, for i = 1, n 
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APPENDIX B DATA SIMULATION 


Simulated pseudorango and delta psoudorango measurements from GPS broad- 
casts are provided by using the ANALYSIS program of GTDS. Further infor- 
mation on the ANALYSIS program is available in Reference 7 . Data may bo 
simulated for either ONPAC or the FILTER and DC programs In GTDS. 

For a data simulation computation, tlie user satellite orbit is computed using 
the GTDS EPHEM program, and the true coordinates are read as input to the 
ANALYSIS program. The GPS satellite coordinates at a given time are deter- 
mined depending on tlie configuration chosen. The options available are 3 GPS 
satellites in 2 orbits (Phase I), G satellites spaced in 3 orbits as set by the 
user (modified Phase I)} 12 satellites Mtli unequal separation (Phase Ila), 

12 satellites with equal separation (Phase Ilb), and 24 satellites with equal sen- 
amtion (Phase III). 

The visibility of the GPS satellites Is constrained by the observer antenna cone 
angle and the Ionospheric cutoff. The user may cycle through all visible GPS 
satellites, select a subset of four by choosing those four which minimize the 
Geometric Dilution of Precision (GDOP), or use all those visible except only 
one of die hvo or more that are essentially In the same direction from the user. 

The user and GPS satellite positions can have random noise added. The GPS 
position and velocity errors are simulated by computing errors in HCL coor- 
dinates, ^vith the along-track error increasing with time. 

The clock errors are a total of the GPS and user satellite clock errors. The 
GPS clock offset Is a constant for each GPS satellite. The constant value Is 
chosen at random for an input standard deviation. The user clock offset may 
be computed using a quadratic or a Markov process model. 
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B 

Tho quadratic model of the user clock offset, , at time t Is 

T®(t) = T® + (t - t^) + 1 £° (t - (B-1) 

B 

where T = tho user clock bias 
u 

B 

f “ the user clock drift which equals tho user frequency offset in parts 
per reference frequency at time t 

f^ = the user clock drift rate which equals tho user discrete change in 
^ frequency per reference frequency 

t « the observation time (in seconds) 

t^ « the epoch of tho clock model 

In the Markov process, the drift rate f^ is computed as a random number at 
discrete update times. The time interval between updates is set by Input, as 
are the mean and standard deviation of the random number generator. The 
drift rate is assumed constant between update times. The clock drift at time t. , 

^ <V’ 


(t.) « (t ) +£ u (t, - 

u 'i' u 'o^ ^ ^ 


(B-2) 


where u(r) is a step function 
(0, T< 0 
T^O 


and the 



are the times of updates. 


The user clock bias is then 


For the quadratic model, the user offset Is computed by evaluating Equa- 
tion (B-1). In ONPAC, the user clock bias, b , Is expressed In kilometers, 

B B 

and Is compared to by multiplying times the speed of light In km/seo. 

The pseudorangc observation is tlie true range plus the user and GPS clock 
offsets plus tlio range measurement noise, computed as a random number witli 
a user-supplied standard deviation! The delta pseudorange Is the difference of 
two true ranges plus Its measurement noise computed from the user-specified 
standard deviation. The filter programs are also supplied with GPS coordinates 
wliich will bo obtained from the broadcasts. The coordinates are Cartesian 
ECEF and contain the simulated GPS ophemeris errors. 

The ONPAC data tape also includes a record with the true observation, true 
user coordinates, and true user clock offset and drift for each observation 
record. This information is used in analyzing the results. 
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Satellite orbit inclination defined in Equation (3-30) 

Imaginary part of the potential expansion defined after 
Equation (3-3) 

Zonal harmonic coefficients 
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Kalman gain vector defined in Equation (2-31) 
Line-of-sight unit vectors defined in Equation (4-3) 
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Satellite mass 
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Covariance matrix defined In Equations (2-33) and 
(2-34) 

Legendre polynomials 

Semilatus rectum defined in Equation (3-28) 

Filter tunable parameter used to compute smoothed 
residual 

State noise process covariance matrix defined in Equa- 
tion (5-G) 

Submatrix components of Q(t) 


Position-velocity and drag submatrix of Q(t) defined 
in Equation (5-8) 

Position-velocity submatrix of Q(t) defined in Equa- 
tion (5-17) 

Component matrices in Q" defined after Equation (5-17) 
Pseudorange defined in Equation (4-1) 

Delta pseudorange, defined in Equation (4-4) 

Mean Earth radius 
Satellite position vector 
Satellite position error 

Real part of potential expansion defined after Equa- 
tion (3-3) 

Satellite height scaling factor, defined after Equation 
(3-6) 

Partial sum in the U-D update defined in Equation (2-30) 
Tesseral harmonic coefficients 
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Gooi’dlnato in Pines force model defined in Equation 
(3-3) 

GPS satellite position components 
Time 

Coordinate in Pinos force model defined In Equation (3-3) 

Filter stepsize 

Maximum allowable stepsize 

Upper right triangular matrix component of the co- 
variance 

Coordinate in Pines force model defined in Equation (3-3) 
Crosstrack unit vector defined in Equation (3-21) 

Element of U matrix 

Gravity potential defined in Equation (3-3) 

Satellite velocity vector 
Satellite velocity error 

Along track unit vector defined in Equation (3-22) 

Diagonal matrix in U-D update derivation 

Diagonal element of W matrix, defined in Equation (2-24) 

Unit vector in UVW coordinate system defined in Equa- 
tion (3-23) 

Vectors of constants used in atmospheric density com- 
putation, defined after Equation (3-6) 

State vector defined in Equation (3-1) 

Error in state vector defined in Equation (5-2) 
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Velocity component of satellite 
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Position component of satellite 
> Velocity component of satellite 
Position component of satellite 
Velocity component of satellite 

Scalar divisor for EKF update defined in Equation (2-15) 
Vectors of smoothed residuals (2-42) 

Filter tunable parameters 

Term in atmospheric density computation defined after 
Equation (3-6) 

Earth's ellipticity 

Atmospheric density defined after Equation (3-6) 

Gravitational parameter, GM 

Fading memory factor, defined in Equation (2-35) 

Ration of the square of the residual to a defined in 
Equation (2-40) 

Tunable parameter used to limit 

Acceleration noise variance rate 
Clock white noise Allan variance factor 
Clock flicker noise Allan variance factor 
Clock random walk Allan variance rate factor 
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Drag noise variance rate 

Measurement variance 

Pseudorange measurement variance 

Delta pseudorange measurement variance 

Time constant of drag defined in Equation (3-12) 

Scaled minimum time where flicker noise predominates 
on the Allan variance curve. 

State transition matrix defined in Equation (3-37) 

Submatrix component of $ defined after Equation (3-37) 

Submatrix component of * defined in Equation (3-39) 

SubmatrLx component of $ defined in Equation (3-38) 

Rotation rate of Earth in Section 3.5.1 

Longitude of the ascending node in Section 3.5.3 

Argument of perigee 
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